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RECENT  PROGRESS  IN  THE  p  AND  h-p  VERSIONS  OF  THE  FINITE  ELEMENT  METHOD 


by 

I .  BabuSka 

Institute  for  Physical  Science  and  Technology 
University  of  Maryland,  College  Park,  MD  20742 


1 .  INTRODUCTION 


-^)The  finite  element  method  has  become  the  main  tool  in  computational 
mechanics.  The  success  is  manifested  by  the  development  of  over  five 

-NA 

hundred  user-oriented  finite  element  program  systems+^see  c?g».  [333^.  The 

literature  on  the  subject  is  overwhelming.  To  date  there  are  over  two 

s’  y 

hundred  monographs  and  conference  proceedings^PWT  and  new  monographs  and 


proceedings  are  continuously  appearing.  Various  forms  of  the  finite  element 
method  are  used  in  practice  for  the  nianerical  treatment  of  elliptic, 
parabolic,  hyperbolic,  linear  and  nonlinear  partial  differential  equations. 


integral  and  integrodifferential  equations,  etc.  Any  class  of  problems  has 
its  own  specific  features.^Xn^this  paper will  only  Steals  with  the  class 

r* 

of  partial  differential  equation  of  elliptic  type.  For  the  sake  of 

iu  — <r^ 

simplicity  we- will  elaborateson  a  characteristic  model  problem  and 


illustrative  results  and  makejonly  additional  comments  of  more  general 


nature  although  the  results  we-fePer  to  are  general. 

hc.fCrfzA-' 


2.  THE  MODEL  PROBLEM  AND  ITS  PROPERTIES 


We  restrict  ourself  to  the  most  simple  model  problem.  Let  (2  c  x  » 

(x^ ,x2)  be  a  bounded  (simple  connected)  domain  with  the  boundary  -  r 
consisting  of  simple  arcs 


flfl(t). 


fif2(t),  t  €  1} 


where  I  -  (-1,1)  and  r  -  U  r, .  We  will  assune  that  f.  ,  and  f,  , 

_i-1  1  1,1  1,2 

are  analytic  functions  on  i  and  fj j  +  f[^2  >  a  >  0.  The  vertices  of 


will  be  denoted  by  Aj ,  i  -  1 . M  and  the  internal  angles  at  Aj  by 


Let  us  be  interested  in  the  model  problem  and  its  standard  (weak 
solution): 


-Au  -  f  on  Q 


(2. 1) 


ropy 

iMtprCTEO  . 


\  \  A  n  >  * 


Here  Q  is  a  subset  of  {1 . M}  -  M.  For  simplicity  we  assume  that 

r°  *  0. 

The  performance  of  any  numerical  method  strongly  depends  on  the 
properties  of  the  (exact)  solution  of  the  solved  problem  especially  on  its 
smoothness.  The  more  information  is  available,  the  better  method  could  be 
desi gned . 

It  is  very  advantageous  to  characterize  the  set  of  solutions  of  (2.1), 
(2.2)  under  consideration  in  the  terms  of  countably  normed  spaces. 

p  p  p 

Let  -  (x1  ifx2  p  and  r£  -  (xj-x^  i )  ♦  (x2~x2  j)  .  Define  *g(x) 

M  B 

-  TT  *V  (x),  B  -  (0fS...9 B|g|) •  0  <  <  1,  and  for  any  integer  k  let 

1-1  M  Bi±k  2 

*6±k(x)  -  r1  (x).  Then  we  let  8g(Q)  «  {u  €  H'(n)| 

l*B+k-2DaulL  (Q)  *  Cdl<k!  •  k  -  2,3.....  |a|  -  k,  C,d  independent  of  k}. 

If  u  €  8o(Q)  then  it  is  analytic  on  Q  -  U  A.  and  has  specific  behavior 

i-1 

in  the  neighborhood  of  Aj,  i  -  1,...,M.  In  [7]  C8]  we  have  proven 

Theorem  2.1.  Let  f  be  analytic  on  Q,  <p  be  analytic  on  r j  ,  J  €  Q  and 

continuous  r°,  ♦  be  analytic  on  r j ,  J  €  M-  Q.  Then  for  (2.1)  and 

(2.2)  u  €  8g(Q)  with  B^  >  B^;  Bj  depends  on  oij  and  the  type  of 
boundary  condition  on  and  Tj.  □ 

If,  for  example,  M  -  5,  Q  -  (1,2,3)  we  get  B^  ■  max(o,1  -  -jj—), 

B2  -  max(o,  1  -  ^J-),  B3  -  max(o,1  -  ,  B4  -  max(o,1  -  B5 

-  max(0,1  -  2J-). 

5 

Remark  2. 1  ♦  In  [9]  we  precisely  characterized  the  traces  of  functions 

p 

from  8|(8)  and  gave  full  characterization  of  the  sets  of  f,  <p  and  41 
which  guarantee  that  the  solution  u  of  (2.1)  (2.2)  belongs  to  8g. 

Remark  2.2.  Theorem  2.1  also  holds  when  the  differential  equation  in  (2.1) 


has  analytic  coefficients  on  S  (see  [7]). 

Remark  2.3.  The  eigenfunctions  of  the  eigenvalue  problem  (2.1)  (2.2)  also 
belong  to  8|(q)  (see  [10]). 

Remark  2.4.  Theorem  2.1  also  holds  for  strongly  elliptic  system  of 
differential  equations  as  elasticity  equations  (see  [11]). 

In  practice,  e.g.  in  the  field  of  structural  mechanics,  the  problems  of 
partial  differential  equations  are  characterized  by  piecewise  analytic  data 
and  hence  theorem  2.1  is  very  well  suited  for  the  applications. 

3.  THE  FINITE  ELEMENT  METHOD 

There  are  various  different  forms  of  the  finite  element  method.  We  will 
consider  here  only  the  basic  class  of  finite  element  methods  (for  our  model 
problem) . 

Let  T  -  {tj}  be  a  partition  of  Q  into  (in  general  curvilinear) 

triangles  or  quadrilaterals  called  elements  t^.  In  the  case  when  T  is  a 

triangulation  we  are  making  the  standard  assumptions.  For  the  general  case 

we  refer  to  [8],  [9],  [36].  We  will  formulate  here  the  assertions  in  the 

case  of  triangulation  only  although  they  hold  in  general.  Let  H(p,T)  - 

{u  €  H1(Q)|u|t^,  €  T  is  a  polynomial  of  degree  p)  be  the  finite 

element  space.  If  tj  is  a  rectangle,  then  polynomials  are  of  degree  p 

in  both  variables.  If  the  elements  are  curvilinear,  then  uL  are  the 

Ti 

standard  "pull-back"  polynomials. 

Remark  3. 1 .  We  have  assumed  that  the  degree  of  polynomials  are  the  same 
over  all  elements.  The  theory  is  developed  for  general  case  when  the 
degree  p  can  be  different  on  different  elements. 

Let  H0(p,T)  -  H(p,T)  fl  Hq(Q)  where  Hq((J)  -  {u  €  H1 ( Q) |u  -  0  on  r°) 
and  Hq( p,T)  be  the  restriction  of  H(p,T)  on  r°.  We  will  assume  that  a 
projection  operator  Pq(p,T)  of  function  tp  into  HQ(p,T)  be  given  and  we 
denote  <pPfj  ■  Pq(p,T)<p.  A  concrete  possible  form  of  Pq(p,T)  will  be 
given  later. 

The  finite  element  method  consists  now  in  finding  upE  -  u(p,T)  € 

H(p,T)  so  that 


1)  u(p,T)  -  <pp>T  on  r° 

2)  //  [pi  |2-  +  |v_)d 

"  '“3x1  3x1  3x2  3x2  1  1 

-  J1  4»v  ds  ♦  jj  fv  dx^  dx2 
r1  Q 

€ 

holds  for  any  v  (  H0(p,T). 

We  will  be  interested  in  the  accuracy  of  the  finite  element  solution 
measured  in  the  energy  norm.  Define  e  -  u  -  Upg  and  let 

i_/-r  _>|2  ff  rrie  \2  .  r?e  \2\.„  . 


|e(T,p)|g  -  JJ  ((|2-)2 .»  (H-)2Kdx 


be  the  error  measured  in  the  energy  norm. 

Two  kind  of  operators  P0(T,p)  can  be  considered.  Let  Y  c  r°  be  the 
side  of  the  triangle  x  with  endpoints  A,B  and  assume  that  Y  ■  [-1,1]. 
A  -  -1,  B  -  1.  Then  <Pp ,7- 1 y  “  where  l(x)  is  linear 

function  on  Y  such  that  <ppj.|Y(±1)  -  (p(±1)  and  a)  in  the  case  of  H1- 

-  P-1 

projection:  Pi(i,p):  <p '  -  I  w'(x)  -  l  a.i’  ,  where  are 

k-0  0  k  k 

the  Legendre  polynomials;  b)  in  the  case  of  H^-projection  P  (T.P);  <p(x) 


l  anVx> 


,  w(x)  -  ^  a, 


kTk(x),  where  Tk(x)  are  the  Tchebyschev 


polynomials. 

As  we  have  seen  there  is  a  large  freedom  in  the  selection  of  H(p,T) 
namely  the  degree  p  and  the  partition  T  and  in  the  selection  of  the 
operators  Pg(T,p).  We  expect  that  |e|E  *0  if  dim  HQ(p,  )  ■*  •. 

It  is  convenient  to  distinguish  in  this  context  three  versions  of  the 
finite  element  method. 

a)  The  h-  vers  ion.  Here  a  sequence  (family)  Hfp,!^)  is  considered 

when  p  is  fixed  (usually  p  -  1,2)  and  the  mesh  Tj_  is  successively 
refined  so  that  the  size  h  of  the  elements  of  goes  to  zero. 

b)  The  p- version.  Here  the  mesh  T  is  kept  fixed  and  p  •*  • 
uniformly  on  selectively. 

c)  The  h-p  version.  In  this  version  the  mesh  is  simultaneously 


refined  and  the  degree  p  increased  uniformly  or  selectively. 

The  h-version  of  the  finite  element  method  is  the  standard  one  and 
extensive  literature  is  devoted  to  it.  The  p-version  is  a  recent 
development.  The  first  theoretical  paper  about  the  p-version  [26]  and  the 
h-p  version  [6]  appeared  in  1981  and  various  results  were  obtained  since 
then.  Some  of  them  will  be  mentioned  later. 

There  are  many  codes,  research  and  commercial  utilizing  the  h-version. 

The  only  commercial  code  using  p  and  h-p  versions  is  the  code  PROBE  which 
was  developed  recently  by  NOETIC  Technologies,  St.  Louis  [5*1].  PROBE  solves 
two  dimensional  problems  of  linear  elasticity,  stationary  heat  problems  and 
thermoelasticity  problems.  The  three  dimensional  extension- of  PROBE  will  be 
released  in  1988.  Three  dimensional  research  code  STRIPE  was  developed  by 
Swedish  Aeronautical  Research  Institute,  see  e.g.  [1].  These  codes  have  in 
addition  various  features  as  adaptive  approaches,  various  a-posteriori  error 
estimation,  etc. 

4.  THE  BASIC  PROPERTIES  OF  THE  p  AND  h-p  VERSIONS  OF  THE  FINITE  ELEMENT 
METHOD 

After  1980  the  p  and  h-p  versions  of  the  finite  element  method  was  studied 
in  detail  from  various  point  of  view.  We  will  mention  here  sane  essential 
illustrative  results. 

In  one  dimensional  setting  the  versions  were  analyzed  in  detail  in 
[34].  Here,  among  other,  the  optimal  meshes  and  p-distribution  has  been 
established  with  upper  and  lower  bounds  of  the  errors  for  the  three  basic 
finite  element  versions. 

In  two  dimensional  setting  the  following  theorem  is  characteristic  for 
the  performance  of  the  h-p  version.  (For  details  see  [8],  [9],  [36].) 

Theorem  4.1.  Let  the  solution  u  of  the  problems  (2.1),  (2.2)  belongs 


to  the  set 

8|(Q).  Then 

there  is  a  sequence  of  meshes 

and  the 

degrees  pt 

such  that 

3 

-o/nT 

M. 

<  C  e 

,  a  >  0 

(4.1) 

where  N^  - 

dim  H0(p1,Ti) 

is  the  number  of  degrees-of- 

freedom  for  the  h-p 

version. 

□ 
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ho tfVOvM 


In  one  dimension  the  rate  is  Ce  It  has  been  proven  in  [34]  that 

the  optimal  mesh  is  a  geometric  one  with  the  factor  (/?  -  1)  -  .17.  The 

experience  shows  that  the  geometric  mesh  with  the  factor  ■  .15  is  also 
optimal  in  two  dimension.  Theorem  *4.1  together  with  Theorem  2.1  shows  that 
practically  in  any  problem  of  structural  mechanics  the  exponential  rate  of 
convergence  can  be  achieved. 

For  the  p-version  the  following  theorem  is  another  typical  one  (for  more, 
see  [21],  [22],  [23],  [24]). 

Theorem  *4.2.  Let  us  consider  the  problem  (2.1)  (2.2)  and  let  B  -  maxd^) 
given  in  Theorem  2.1.  Then  for  the  p-version  we  have 

|e|E  <  (*4.2) 

while  for  the  h-version  with  uniformed  mesh 

-(III) 

|e|E  >  CN  2  .  (4.3) 

□ 

In  Theorems  4.1  and  4.2  either  the  projection  Pq  or  P^  could  be  used 
provided  <p,  i>  are  sufficiently  smooth.  The  difference  between  these  two 
operators  occurs  for  the  p-version  when  the  boundary  condition  <P  is 
unsmooth,  e.g.  cp  €  H*(r),  V2  <  6  <  3/4.  In  this  case  the  projection 

has  to  be  used.  For  the  h-p  version  there  is  no  difference  in  the 
asymptotic  rate  but  some  difference  occurs  in  the  constant  of  the  estimates. 
For  the  analysis  of  the  influence  of  the  operator  PQ  on  the  accuracy  we 
refer  to  [12]. 

Remark  4. 1 .  So  far  we  have  assumed  that  the  domain  (2  is  bounded. 

Nevertheless  the  exponential  rate  of  convergence  (4.1)  holds  also  for  the 
r  2 

problem  on  Q  -  R  -  Q  when  f  has  bounded  support.  Here  the  infinite 
elements  and  properly  selected  shape  functions  have  to  be  used.  For  more, 
see  [15]. 

Remark  4.2.  The  h,  p  and  h-p  versions  have  different  aspects  with  respect 
to  the  pollution  problem.  In  presence  of  a  singular  behaviour  of  the 
solution  (e.g.  in  the  neighborhood  of  the  entrant  corner  of  the  domains) 


the  L,,,  error  is  very  large  in  the  element  consisting  the  corner.  This 
effect  disappears  in  elements  which  are  separated  away  from  the  singularity 
by  few  elements.  This  effect  is  essential  for  a  proper  mesh  design  in 
practical  computation.  For  details  we  refer  to  [14]. 

So  far  we  have  dealt  with  the  problem  (2.1)  and  (2.2)  of  second  order. 
For  the  analyses  of  the  finite  element  solution  for  the  problems  of  order 
2k  we  refer  in  the  case  of  the  h-p  version  to  [35]  and  the  p-versi  n  to 
[39]  [53].  For  the  basic  analysis  of  the  p-version  in  3  dimensions,  we 
refer  to  [30]  and  [31].  The  eigenvalue  problem  is,  as  is  well  known, 
directly  related  to  the  ‘•source"  problem  we  addressed  earlier.  See  e.g. 
[16]  and  [17].  In  the  case  of  the  eigenvalue  problem  (in  our  case) 

Au  -  Xu 

u  -  0  on  r° 


fa  -  o 

9u 


on 


,(1> 


The  eigenfunctions  belong  to 


|  x-x 


FE I 


S  C  e 


3 

-2a/N 


5g(n)  and  hence 


luFE^Xh^  ”  u^Ie  -  C  e 


3 

-a/N 


For  more  details,  see  e.g.  [10]. 


5.  IMPLEMENTATION  AND  COMPUTATIONAL  COMPLEXITY 

There  are  some  essential  features  of  implementation  of  the  p  and  h-p 
versions.  See  e.g.  [56],  The  elements  are  of  a  hierarchical  type  which 
leads  to  augmentation  (bordering)  of  the  local  stiffness  matrices  when  p  is 
increased.  This  also  allows  to  change  very  flexibly  the  degree  of  the  shape 
functions  from  one  element  to  the  other  one.  The  shape  elements  are  (in  two 
dimensions)  of  nodal  type,  side  type  and  internal  type  and  are  based  (in 
PROBE)  on  the  integrals  fo  the  Legendre  polynomials.  This  is  important  for 
numerical  stability  aspects.  The  computation  (in  two  dimension)  of  the 
local  stiffness  matrix  on  a  rectangle  requires  (when  optimally  programmed) 
0(p4)  operations.  The  computation  of  the  local  stiffness  matrices  takes 
much  more  effort  for  the  p-version  (with  high  p)  that  for  the  h-version  and 


hence  the  p-version  is  very  well  suited  for  parallel  computations.  The 
sparsity  of  the  global  stiffness  matrix  is  also  smaller  for  the  p  version 
than  the  h-version.  Hence,  the  complexity  of  the  computation  for  the  same 
number  of  degrees-of-freedom  is  higher  for  the  p-version  than  the  In¬ 
version.  Nevertheless,  it  is  essential  to  relate  the  achieved  accuracy  to 
the  computational  effort.  For  an  analysis  of  the  computational  complexity 
and  computer  time  comparison,  we  refer  to  [19]  and  [20].  The  results  show 
that  the  h-p  version  with  higher  degree  p  is  preferable  for  solutions 
which  are  not  overly  unsmooth  or  have  singular  behaviour  in  a-priori  known 
areas  as  in  the  neighborhood  of  the  corners.  If  the  solution  is  very 
unsmooth,  e.g.  if  the  coefficients  in  the  differential  equations  are 
uniformly  rough,  for  example,  measurable  only,  then  only  low  accuracy  i3 
practically  achievable  by  any  method  and  h-version  with  low  p  is 
preferable.  In  general,  for  very  low  accuracies  the  low  order  elements  are 
preferable,  for  the  modest  one  higher  degree  elements  have  to  be  preferred. 

6.  THE  PROBLEM  OF  THE  MESH  DESIGN 

One  of  the  most  laborious  part  of  the  finite  element  analysis,  especially  in 
three  dimensions,  is  the  mesh  generation  also  when  sophisticated  mesh 
generators  are  used.  The  use  of  large  elements  (possibly  of  high  degrees) 
which  are  describing  only  the  geometry,  greatly  simplifies  the  user's  work, 
also  if  possibly  on  expense  of  the  computer  time.  (It  Is  necessary  to 
realize  that  the  relation  between  manpower  cost  and  computation  cost  is 
going  steadily  up.)  The  option  of  a  change  of  the  degrees  of  elements 
Increases  significantly  the  flexibility  of  the  program  and  gives  the  user 
effective  tool  for  the  quality  control.  The  p  and  h-p  version  programs  give 
such  options.  It  is  advantageous  therefore  to  create  directly  or  indirectly 
the  proper  mesh  and  to  achieve  then  the  desired  accuracy  by  an  appropriate 
choice  of  the  degrees,  which  can  be  made,  for  example,  in  an  adaptive 
mode.  The  goal  is  to  achieve  the  same  combination  of  the  degrees  and  mesh 
refinement  which  would  be  obtained  for  the  given  acctracy  by  the  h-p  version 
directly.  To  achieve  this  goal  two  avenues  could  be  followed,  the  expert 
system  and  the  adaptive  approach.  The  expert  system,  see  e.g.  [18],  [45], 
advises  the  user  how  to  design  the  mesh  and  element  degrees  for  the 
requested  accuracy  and  provides  the  user  with  a  mesh  generator.  The  expert 
system  is  interactive,  follows  the  progress  of  computation,  gives  the  user 


on  his  request  various  desired  information  for  an  effective  computation  and 
engineering  analysis.  The  experience  (see  [18]  [19])  is  that  the  cost  of 
the  expert  activity  for  a  mesh  and  degree  design  is  at  most  20%  of  the  total 
cost. 

The  adaptive  approach  (see  e.g.  [46],  [34],  Part  3)  which  is  possible  to 
see  as  an  "automatic"  expert  system  makes  various  decisions  for  the  users. 
Both  approaches  have  some  common  parts  but  the  concepts  are  significantly 
different.  We  refer  also  to  [46]  for  various  additional  aspects. 

7.  THE  ROBUSTNESS 

An  effective  method  has  to  perform  uniformly  well  for  a  broad  class  of  input 
data.  The  elasticity  problems  can  be  in  practice  nearly  singular  as,  for 
example,  in  the  case  of  nearly  incompressible  material,  various  plate  and 
shell  theories,  in  the  case  of  thin  domains,  etc.  The  h-version  suffers  in 
these  cases  by  the  "locking"  problem  which  has  to  be  overcome  by  various 
special  approaches  as  reduced  Integration,  etc.  (see  e.g.  [28]).  Problems 
of  these  types  are  usually  avoided  when  the  p  and  h-p  versions  are  used. 

The  convergence  rate  then  (ir.  contrast  to  the  h-versions)  is  uniform  with 
respect  to  the  Poisson  ratio  when  higher  degrees  of  elements  are  used.  See 
e.g.  [49],  [50]  [64]. 

8.  THE  QUALITY  CONTROL  OF  THE  SOLUTION 

It  is  essential  to  have  a  possibility  of  a  quantitative  assessment  of  the 
quality  of  computed  data.  For  the  survey  of  today's  general  ideas  and 
results  in  this  direction,  we  refer  to  [41].  In  the  case  of  the  p  and  h-p 
versions  there  is  relatively  easy  way  for  the  quality  assessment  of  any  data 
of  the  interest  by  changing  the  degrees  and  by  an  extrapolation  procedure. 
This  approach  is  very  effective  because  it  indicates  reliably  the  errors  of 
any  computed  data  of  interest,  the  energy  norm,  value  of  stresses,  stress 
intensity  factor,  etc.  See  e.g.  [1],  [8],  [36],  [58],  [61]. 

9.  THE  COMPUTATIONAL  AND  ENGINEERING  EXPERIENCE 

Because  of  the  developed  commercial  code  PROBE  and  research  code  STRIPE  an 
extensive  experience  is  already  available  in  the  research  and  industrial 
use.  For  the  industrial  experience  we  refer  e.g.  to  [27]  where  numerical 
results  are  presented.  For  the  research  one  we  refer  e.g.  to  [1],  [3],  [8], 


C9],  CIO],  [18],  [19],  [20],  [25],  [36],  [48],  [55],  [57],  [63].  The 
experience  shows  that  the  p  and  h-p  versions  of  the  finite  element  method 
has  many  practical  advantages  for  the  engineering  computations  for  linear 
elliptic  problems. 


10.  RELATION  TO  SOME  OTHER  METHODS 

The  ideas  of  the  p  and  h-p  versions  are  related  for  example  to  the  various 
methods  used  in  fluid  dynamics  as  the  spectral  method  and  its  variation.  We 
refer  e.g.  to  [43]  and  references  given  there.  A  commercial  fluid  dynamics 
code  written  by  T.  Patera  is  closely  related  to  the  idea  of  h-p  version  of 
the  finite  element  method. 

The  h-p  version  of  course  generates  various  finite  element  approaches  and 
in  principle  encompasses  various  diverse  approaches,  see  e.g.  [29].  The  h-p 
version  can  be~naturally  also  used  for  solving  integral  equations,  boundary 
element  method,  etc.  See  e.g.  [47],  [51],  [67]. 

11.  SUMMARY 

The  p  and  h-p  version  of  the  finite  element  method  is  a  new  development 
which  gives  new  possibilities,  theoretical  and  practical  for  the  finite 
element  method.  It  Is  today  reasonably  well  understood  in  the  case  of 
elliptic  equations,  both  theoretically  and  practically.  The  aim  of  the 
present  paper  was  to  give  a  brief  survey  of  various  aspects  of  the  p  and  h-p 
versions  of  the  finite  element  method  for  solving  elliptic  linear  problems 
and  provide  the  comprehensive  references.  Nevertheless,  many  theoretical 
and  practical  aspects  of  the  method  for  other  problems,  linear  and  nonlinear 
are  still  to  be  resolved  as  well  as  the  problems  of  implementation  for  three 
dimensional  problems. 

In  addition  we  refer  to  the  references  [2],  [3],  [4],  [37],  [38],  [42], 
[44],  [52],  [58],  [60],  [62],  [65],  [66]  which  are  directly  related  to  the 
subject  discussed. 
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